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Abstract 

We study the damage spreading transition in a generic one-dimensional 
stochastic cellular automata with two inputs (Domany-Kinzel model) Us- 
ing an original formalism for the description of the microscopic dynam- 
ics of the model, we are able to show analitically that the evolution of 
the damage between two systems driven by the same noise has the same 
structure of a directed percolation problem. By means of a mean field ap- 
proximation, we map the density phase transition into the damage phase 
transition, obtaining a reliable phase diagram. We extend this analysis 
to all symmetric cellular automata with two inputs, including the Ising 
model with heath-bath dynamics. 

Key Words: Damage spreading, directed percolation, stochastic cel- 
lular automata, disordered systems, symmetry breaking. 
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1 Introduction 

In this paper we deal with the problem of the evolution of two replicas of a 
Boolean system (cellular automaton) that evolve stochastically under the 
same realization of the noise. The system is defined on a regular lattice of 
L sites and evolves in discrete time steps. We limit the explicit analysis 
to one dimensional systems, but the results can be extended to higher 
dimensions. 

Let us indicate the time with the index t = 1, . . . , oo and the space 
with i = 0, 1, . . . , L — 1. The state variables a(i, t) can assume the values 
or 1 (Boolean variables). The evolution of cr{i,t) is given by probabilistic 
transition rules and depends on a small number of neigboring sites; in 
its simpler form, a(i,t) depends only on the state of the two nearest 

* Dipartimento di Matematica Applicata, Universita di Firenze, via S. Marta, 3 1-50139, 
Firenze, Italy; INFN and INFM sez. di Firenze. e-mail:bagnoli@ing. unifi.it 



1 



t(0,0- 


-1) 


= Po 


r(0,l- 




= Pi 


t(1,0- 


-1) 


= P2 


t(1,1- 


-1) 


= £>3 



neigbors. In this case one can consider the space-time lattice as a tilted 
square lattice. 

In order to simplify the notation, we write <r+ = a(i + l,t), cr_ = 
cr{i— l,t), u' = a(i,t+l). The evolution rule can be synthetically written 
as 

a = f(a-,a+). 

Since the number of possible values of the couple (<r_, a+) is four, the 
function / is usually specified by giving the four transition probabilities 
t(<t_,(t+ —* 1) from each possible configuration to one: 



(1) 



The normalization condition gives r(a-,cr+ — ► 0) = 1 — t(ct_,ct+ — > 1). 

All sites of the lattice are generally updated synchronously. Except 
for the case of deterministic cellular automata, for which the transition 
probabilities are either zero or one, we do not expect strong differences 
between parallel and sequential updating. 

This schematization naturally arise in the modelization of several sys- 
tems (contact processes), in physical and biological investigations. It has 
been introduced by E. Domany and W. Kinzel [jj], and can be consid- 
ered the prototype for all local stochastic processes. For a short review of 
the applicability of this model, see references || and |^|. 

In the thermodynamical limit, the Domany-Kinze 
hibits a phase transition from an ordered to a disordered phase for po — 0. 
The ordered configuration is a(i) = for all i (adsorbing state). The order 
parameter is the asymptotic density m = lim t ^oo limL^oo m(t,L), where 
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In the following, we refer to the critical surface m — and all its inter- 
sections with planes in the pj space with the symbol a (see the figures). 

This transition has been studied mainly for the symmetric case pi = 
P2- Except for a phenomenological renormalization study to which we 
collaborated [Q, the transition line has been found numerically to belong 
to the universality class of directed percolation, which is a particular case 
of the model. The disagreement for the renormalization group results 
can originate from finite-size effects. For the asymmetric case p\ ^ P2, 
it has been claimed ^) that the phase transition belongs to a different 
universality class (mean field). 

The existence of an adsorbing state is a non-equilibrium feature of 
the model, allowing the presence of a phase transition also in a one- 
dimensional (spatial) system. It is shown is section^ that in the DK 
model there can be two adsorbing states (ff(i) = and a(i) = 1), related 
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by a simple transformation of the transition probabilities. The two transi- 
tion lines met at the point M (pi = 1/2, pa = 1). This point corresponds 
to the problem of a random walk in one dimension, and thus exhibit mean 
field exponents. 

A powerful tool for the investigation of this kind of models is the 
study of damage spreading. One considers two replicas a and 77 of the 
same model, with different initial conditions (they can be completely un- 
corrected or differ only in some sites). The two replicas evolve under the 
same realization of the stochasticity. The difference at site i and at time 
t between the two configurations is given by 



where the symbol © represents the sum modulus two (eXclusive OR, 
XOR). Since we use Boolean variables (a, b G {0,1}) one can interpret 
the exclusive oiasa(Bb — a + b — 2ab. When mixing XOR and AND 
(represented as a multiplication), one can use the algebraic rules for the 
sum and the multiplication. 

The order parameter for the damage spreading transition is the asymp- 
totic Hamming distance H = lim t _, 00 linii^oo H(t, L) defined as 



using the usual sum. 

The critical surface H = and its intersection are indicated with the 
symbol 7. 

In the DK model, numerical and analytical investigations j(| ^, |^, 
[lo| , indicated the existence of a damage spreading phase. 

The damage phase transition can be thought as an ergodicity breaking 
transition: in the phase where the damage disappears, all initial conditions 
asymptotically follow a trajectory that does not depend on the initial con- 
ditions, but only on the realization of the noise. The Hamming distance 
can be easily related to the overlap between the configurations. 

The critical exponents for the density and the damage transitions in 
the planepi = P2 and po ~ are numerically the same [| . It has been 
guessed pj] , [l2"[ that all continuous transition from an adsorbing to and 
active state belong to the universality class of the DK model (and thus of 
directed percolation), and that the same universality class should include 
all damage spreading transitions j^] . 

Here we want to investigate the connection between the density phase 
transition and the damage phase transition in the DK model. We have to 
carefully describe the dynamics of the model: the position of the transition 
line depends on the way in which the randomness is implemented in the 
actual simulations. In section ^| we introduce the formalism that allows 
an exact description of how randomness is implemented in the model. We 
are thus able to write down the evolution equation for the spins, and to 
obtain the evolution equation for the distance between two replicas. The 



h(i,t) = a{i,t)®r]{i,t); 
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structure of the latter equation corresponds to the DK model with po — 0. 
We conclude that the universality class of damage spreading is, at least 
for this simple case, that of directed percolation. In section ^ we obtain 
the phase diagram of the DK model by mapping the transition line for 
the density to the transition line for the damage by means of mean field 
approximations. In section ^| we show that one can infer the existence 
of a phase transition for the damage also in cases for which there is no 
phase transition for the density, and that there are two disjoint regions in 
the parameters space for the damage spreading. Finally, conclusions and 
open questions are drawn in the last section. 

2 The damage spreading transition 

Let us start from a simple example, the dilution of rule 90 (in Wolfram's 
notation that will also serve to fix the notation. Rule 90 is a deter- 
ministic rule that evolves according with 

a — a- ® a + . 

The transition probabilities for the diluted rule 90 are 
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p being the control parameter of the model. 

In order to apply rule 90 for a fraction p of sites, and rule (all 
configurations give 0) for the rest, one usually extracts a random number 
r = r(i, t) for each site and at each time step and chooses the application 
of rule 90 or rule according with r < p or r > p resp. 

We can easily write the explicit expression for this rule by means of 
the function [•] , assuming that [logical proposition] takes the value 1 if 
logical proposition is true, and otherwise (this interpretation of logical 
propositions is the standard one in C language). Finally, we have for the 
diluted rule 90 

a = [r < p] {a- 0ct+). (2) 

One can also think of having all random numbers r(i, t) extracted 
before the simulation and attached to the sites of the space-time lattice 
even though they are not always used. The random numbers are thus 
similar to a space-time quenched (disordered) field. 

Once given the set of random numbers, the evolution is completely 
deterministic, and the evolution function depends on the lattice position 
(spatial and temporal) via the random numbers r(i,t). One can alterna- 
tively define the model stating that deterministic functions are randomly 
distributed on the space-time lattice according to a certain probability 
distribution. This description is very reminiscent of the Kauffman model 

& 
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The damage spreading can be considered a measure of the stability of 
the set of possible trajectories, averaging over the realizations of the noise. 
The original definition of Lyapunov exponent is a measure of the instan- 
taneous effects of a vanishing perturbation. Since the state variables of 
cellular automata assume only integer values, one has to extend the defini- 
tion to a finite initial distance (and to finite time steps) , thus taking into 
account the possibility of non-linear effects. For cellular automata, the 
smallest initial perturbation corresponds to a difference of only one site 
between the two replicas. The short-time effects of a (vanishing) pertur- 
bation define the analogous of the derivatives for a continuous system Jl^] . 
The study of the equivalent of the usual (linear) Lyapunov exponent for 
deterministic cellular automata allows a classification of the rules accord- 
ing with the trend of the damage E9|. The general problem of damage 
spreading can thus be considered equivalent to the study of the non-linear 
Lyapunov exponent (i.e. finite initial distance and finite evolution times) 
for space-time disordered cellular automata. 

Using the concept of Boolean derivatives j^], we develop a Boolean 
function f(a,b) as 

f(a,b) = f ® fia® frb® fcab, 

where the Taylor coefficients are 

fa = /(0,0) 

fi = /(0, 1)0/(0,0) 

h = /(l, 0)9/(0,0) 

h = ,f(l, 1)9/(0, 1)©/(1, 0)©/(0,0). 

One can verify the previous expression by enumerating all the possible 
values of a and b. 

Using the bracket [•] notation, the transition probabilities (|l|) corre- 
spond to 

/(0,0) = [r <p ] 
/(0, 1) = [ri < Pi] 
/(1,0) = [r 2 < P2 ] 
/(1,1) = [r 3 <p 3 ], 

where the random numbers rj(i,t) belongs to the interval [0, 1) and con- 
stitute the quenched random field. We neglect to indicate the spatial and 
temporal indices for simplicity. 
The Taylor coefficients become 

fo = [ro < Po] 

fi = [n < Pi] © [r < po] 

/a = I r 2 < Pa] ® [ro < Po] 

/a = [r 3 < ps] © [r a < Pa] ffi [n < pi] ® [r < Po]- 

In the following we shall assume pi = p 2 and n = r 2 so that /i = f 2 
and 
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h = [r 3 < p 3 ] © [ro < Po]- 

The correlations among the random numbers Tj (at same space-time 
position) affect the position of the critical line for the damage (V). as 
already pointed out also by P. Grassberger || and E. Domany |l7l] , but 
not the position of the transition line for the density (a). Only a careful 
description of how the randomness is implemented in the model completely 
specify the problem of damage spreading. In principle one could study 
the case of generic correlations among these random numbers. Here we 
consider only two cases: either all rj are independent (case i, transition 
line 7i) or they are all identical (case m, transition line yu). 

The evolution equation for the single site variable a — a(i, t) is 

a' = [r < po] © ([n < pi] © [ro < po]) {a- © 0+) © 

([r 3 <p 3 ]©[r <p ])ff_ff+. (3) 

We can substitute the evolution equation for the replica 7} = a © h, 
with the evolution equation for the damage h — a © 77, obtaining 

h' = ([ri <pi]ffi[r <po]®([r 3 < p 3 ] © [r < p ]) er+) h_ © 
([ri < pi] © [r < po] © ([r 3 < Ps] © [ro < Po]) cr_) fe+ © 
( [r 3 < pa] © [ro < Po] ) ft- h+ . (4) 

This equation has the same structure of the evolution equation of the 
original model Q if in this latter we set po — 0. Remembering that only 
for this value of po the DK model exhibits a phase transition, we have 
a strong argument for the correspondence between directed percolation 
and damage spreading transitions. However, also in the symmetric case 
Pi = P2 and ri = V2, the evolution equation of h is symmetric only in 
average, and one has to take into consideration the correlations between 
cr_ and <r+ . As discussed before, these correlations can be included in the 
definition of the DK model, which specify only the transition probabilities. 
It remains to be proved that all these versions do belong to the same 
universality class. 

Let us now consider the case po = 0. Previous numerical investigations 
showed that on this plane the two curves a and 7 meet at the point 
Q = (~ 0.81, 0). Inserting the value p 3 = po = in the equation (E), we 
see that the evolution law for h is the same of that for a, and so both 
transitions coincide on this line. This corresponds also to the dilution of 
rule 90. 

Since the rest of the 7 curve lies away from the density transition line, 
the correlations among sites decay rapidly in time and space. This allow 
us to use a mean field approximation. In the simplest form, we replace 
a(i,t) with a random bit that assumes the value one with probability m. 
With this assumption the equation (Q) becomes 

ti = ([n < pi] © [7-3 < p 3 ][r 4 < m])h- © 

([ri < pi] © [r 3 < p 3 ][r 5 < m])h+ © [r 3 < p 3 ]h-h+; (5) 
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where m and rs are independent random numbers. This is a strong ap- 
proximation, both because of correlations and because the same o(i,t) is 
shared by h(i — l,t + 1) and h(i + 1,4 + 1). Nevertheless, we can assume 
this equation as a starting point in our derivation of the phase diagram. 

We now want to remap this model onto the original DK model, as- 
suming that the asymmetry {n 7^ rs), that in average vanishes, does not 
strongly affect the transition. 

The remapped transition probabilities p are 

f(0,0->-l) = po = 

f (0, 1 — > 1) = pi = 7T ([n < pi] © [r 3 < pa] [rs < m}) 

f(l,0— > 1) = pi = n ([n < pi] ® [r 3 < p 3 ][r 4 < m]) 

f(l,l-t>l) = pa = n ([r s < p a ]([r 4 < m] © [r 5 < m] © 1)) , 

where 7r(/(r)) = J Q drf(r) is the probability that the Boolean function / 
of the random number r takes the value one. 
For case i (n =£rz), we have 



pi = pi + pzm - 2p 1 p 3 m 
P3 = ps(l-2m(l-m)), 



(6) 



while for case ii (n = rs) 



(7) 



pi = m|pi -P3I + (1 - m)pi\ 
P'i = ps(l — 2m(l — m)). 

Since 7 lies in the pi > P3 region, one has for case it 
pi = Pi — mp3- 

Notice that for P3 = or for m = the two curves "y» an d 7« coincide, 
as already noticed numerically by Grassberger 

Given a certain point (pi,P3), it belongs to the damage transition line 
7 (H(pi,ps) — 0) if the point (pi,P3) belongs to the density transition line 
a (m(pi,p3) = 0). In order to draw the phase diagram for the Hamming 
distance, one has to know the value of the density m in all the parameter 
space, and in particular the position of a. Unfortunately, we do not have 
a simple expression for these quantities; in the next section we use some 
approximation in order to draw a rough phase diagram. However, we are 
able to demonstrate that a and 7 are tangent at point Q. 

The slope q of the normal to a at Q can be given as 



dm j dm 
dpi / dps, 



Q 



Considering that 7 = H(p!,p 3 ) = = m(pi(pi,p3),p 3 (pi,p 3 )) = 0, 
the partial derivatives of H are given by 

dH dm dpi dm dps 

dpi dpi dpi dp3 dpi ' 

dH dm dpi dm dpz 

dp3 dpi dp3 dp3 dp3 ' 
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One has to take into account that pj depends on pi both directly and 
via m. Inserting the relations ^| or ^ and considering that at point Q, 
m = P3 = one obtains that 



q = 



dH 
dpi , 



dH 

dp 3 



= <?■ 



Since we know from numerical experiments and from all mean field ap- 
proximation beyond the very first one that the slope of a at Q is negative 
in the (jJi,f>3) plane, the tangency of 7 to a implies a reentrant behavior 
for the damage transition curve, as observed in reference tel. 



3 The phase diagram 

The problem of sketching an approximate phase diagram for the damage 
in an analytical way has been dealed with by several authors |jj 
Since any equation for the damage depends on the behavior of one replica, 
there are two sources of errors to be controlled: the approximations for 
the evolution of one replica and that for the difference (or for the other 
replica). As a consequence, all approximation schemes proposed so far 
require large efforts for a poor result. Our method is able to exploit the 
knowledge of the density phase to study the damage phase transition. 
There are several methods that rapidly converge to a good approximation 
of a; to our knowledge the best ones are the phenomenological renor- 
malization group || and the cluster approximation (local structure) (jl8j 
improved by finite-size scaling. This latter method can also give a good 
approximation of the behavior of ra(po,Pi,P3) at any point. 

Since here we are not interested in numerical competitions, we use 
the high-quality data for the density transition line from reference Q 
combined with a first order mean field approximation for the density. 
The a curve has been approximated in the (jJi,P3) plane by a 5 th order 
polynomial 

5 

P3 = }^a,ip\. (8) 

i=0 

The simplest mean field approximation for the asymptotic density m gives 

l-2 Pl 

m = . 

P3 — 2pi 

By using these approximations one obtains from equations (^) or (Q) 
the curves reported in figure [j], together with the presently best numerical 
results ||.The main source of error is that in the mean field approxima- 
tion the a curve does not corresponds to the zero of the density. This is 
particularly evident in the absence of reentrancy of curves 71 and 72. Nev- 
ertheless even this rough approximation is able to reproduce qualitatively 
the phase diagram and to exhibit the influence on the damage critical line 
of the different implementations of randomness. Notice that the 7 curve 
from reference M corresponds to the implementation of equation (pf). 
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4 The po > case 



The DK model with arbitrary p includes all one dimensional symmetric 
cellular automaton model or spin system with two inputs. We can rep- 
resent each possible model as a point in the three dimensional unit cube 
parametrized by po, Pi, P3- The general form of the transition probabili- 
ties from equation ([|) is 

Pi = pi + (l-2pi)(mp 3 + (l-m)p ); / g v 
P3 = (po +P3 - 2p p 3 )(l - 2m + 2m 2 ). 

There is a trivial transformation of the original DK model with po — 0. 
One can revert (0 <-» 1) all the spins before and after the application of 
the rule. The new transition probabilities p\ are 

Po = 1-P35 

pi = i-pi; 

P3 = 1 - PO- 

The critical plane po = maps to P3 = 1, and the adsorbing state is 
now the configuration in which all spins are one. We indicate with the 
symbol a' the critical curve obtained by this transformation. The point 
Q is mapped to the point Q' — (1, ~ 0.2, 1). The parameter cube and 
the critical curves are reported in figure m. This mapping suggests the 
presence of a damaged zone near the corner (1, 0, 1). 

In order to study the position of the critical surfaces for the damage, 
we numerically solved equation Q combined with the expression (^) of 
the critical line a in the very simple approximation m = 0.5. The results 
are reported in figure ^. Direct numerical simulations qualitatively agree 
with this picture. 

The one dimensional Ising model in zero field with heath bath dynam- 
ics can also be expressed with this formalism. 

The local field g = g% for the one dimensional Ising model is 

S = J£-((2<7_-l) + (2<7 + -l)), 

where K = (5 J — J/ksT is the rescaled coupling constant and a = 0, 1 
the site variables (spin). The local field g can assume the values ~2K, 0, 
IK. 

For the heath bath dynamics, a' takes the value one with probability 
p given by 

1 

P ~ 1 + exp(-2 ff ) ' 
The transition probabilities are 

Po = TTI ; 
i 

pi = 5; 
P3 = ITT' 
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where £ = exp(— 4K). Notice that pz = 1 — po; for T > 0, po < 1/2, while 
for negative temperatures po > 1/2. The point po = Pz — 1/2 corresponds 
to infinite T. 

The evolution equation for the site variable is 

a = [r <p ] © ([r <pi] © [r <p ])(a- fficr+) © 
([r < P3] © [r < po])cr_cr + ; 

where usually all Taylor coefficients depend on the same random number 
r — r(i, t). The existence line ui+ for the Ising model with T > 0, p\ = 1/2, 
pz — 1 — Po, intersects the critical line a at A/ = (0, 1/2, 1). The existence 
line to- for T < ends at M' = (1 1/2, 0). The point R = (1/2, 1/2, 1/2) 
corresponds to T — oo (see figure 

The evolution equation for the Hamming distance h is equivalent to 
equation (Q), with all rj equal to r. Taking into account the correlations 
induced by the random numbers, and that the magnetization is 1/2 except 
at the critical point, one obtains 



Pi = 
P3 = 



2(1+0' 

i+r 



i.e. the line x = Pi = 2p2, po = for positive or negative temperatures. 
The line \ intersects a at point M for T = , confirming that the sym- 
metry breaking transition for the Ising model occurs at zero temperature. 



5 Conclusions and perspectives 

In this work we presented a formalism that allows the careful description of 
Boolean algorithms for stochastic cellular automata (including spin system 
like the Ising model) . Using this formalism we were able to derive the exact 
equation for the evolution of a damage between two replicas that evolve 
under the same realization of the noise. Using a mean field hypothesis, 
we gave strong indications that the critical line for the damage phase 
transition belongs to the same universality class of that for the density in 
the DK model, and thus to the directed percolation universality class. We 
mapped the density critical line to the damage critical line, obtaining the 
regions in the parameter space of a general symmetric cellular automaton 
where the replica symmetry breaking is to be expected. Our predictions 
are qualitatively confirmed by numerical simulations. 

Several questions remain to be answered. Among others: is it possible 
to obtain similar results starting from a field description? How does the 
phase diagram for more general (asymmetric, three-input, etc.) cellular 
automata look like? 
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Figure Captions 



1. Phase diagram for the density and the damage in the DK model 
(po = 0). The curve labeled a is the density transition line and the 
one labeled 7 is the damage transition line from reference the 
curves labeled 7* and 7« correspond to mean field approximation of 
equations (|^) and (0) resp. 

2. The parameter cube for the general symmetric cellular automata. 
The dashed curves labeled a and a' belong to planes po = and 
P3 = 1 resp., and correspond to the density phase transitions. The 
solid curves correspond to the intersection of the damage critical 
surface (shaded) 7 and 7' with the boundaries of the cube. The 
dotted-dashed lines labeled lu+ and u- correspond to the existence 
line for the Ising model for positive and negative temperatures, resp. 
The points labeled M and M' to the critical points of Ising model at 
zero temperature, and the point labeled R to infinite temperature. 
The dotted line labeled \ corresponds to the damage in the Ising 
model. 
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F. Bagnoli - damage spreading transitions 




figure 1 



F. Bagnoli - damage spreading transitions 
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